In seismic data acquisition, the seismic source is typically positioned at a selected shot location, and the seismic reflections of the shot are detected (the “shot record”) by receivers also located at selected locations. Then, the source and receivers are moved to different locations and the process repeated, and in this manner a seismic survey is taken of a selected subterranean region. Ideally, for many seismic processing and interpretation objectives, the source and receiver locations would lie on a uniformly and densely sampled grid, but this is difficult to achieve in common industry practice for many reasons including surface obstructions, currents, cable feathering, and acquisition cost. Survey economics mandate that the spacings be as large as will still yield the required detail in the survey results. The desired seismic reflections are wavefields that reflect once from a subterranean interface between regions with different acoustic properties (such as the upper surface of a petroleum reservoir), and then travel back to the surface to be detected by a receiver. This desired data are often obscured by undesired multiple reflections of seismic rays. Multiple suppression techniques exist for reducing this noise problem. In Surface-Related Multiple Elimination (SRME), multiple reflection that have at least one downward reflection from the free surface (interface between air and water or air and land) are predicted (see for example Dragoset, and Jericevic, “Some remarks on surface multiple attenuation,” Geophysics 63, 772-789 (1998)). In contrast, Interbed (Internal) Multiple Elimination (IME) handles multiple reflections with downward reflections in the subsurface (see for example Jakubowicz, “Wave equation prediction and removal of interbed multiples,” 68th Annual International Meeting, SEG, Expanded Abstracts, 1527-1530 (1998)). Some implementations of the latter algorithm require identification of strong reflectors that act as the main interbed multiple generators, which can be accomplished through an application of a VSP deconvolution technique described in Ross and Shah, “Vertical seismic profile reflectivity—Ups over downs (short note),” Geophysics 52, 1149-1154 (1987). Multiple predictions produced by SRME and IME are typically not accurate enough to be directly subtracted from the input multiple-contaminated data. Thus, an additional adaptive subtraction step is required, where predicted multiples are first shaped to better fit the actual multiples present in field data and then subtracted from the input data. Alternatively, multiples in field data can be attenuated based on their closeness to the shaped predicted multiples.
Data-driven multiple prediction methods (SRME and IME) require pair-wise convolution of densely and regularly sampled shot gathers centered at source and receiver positions for each trace. This requirement is the key contributor to the high computational cost of these methods. A straightforward implementation of SRME is illustrated in FIG. 1A and requires the following main steps:                1. Densely sampled shot gathers are reconstructed at each surface location where either a source or a receiver was present at any time during field acquisition. (The term “shot gather” describes a gather with an impulsive pressure source at the center of the gather with pressure sensors (receivers) all around it. Due to seismic reciprocity (i.e., the fact that recorded data would not change if we interchanged source and receiver positions), a “shot gather” can be also thought of as a “receiver gather” with a pressure sensor in the middle and impulsive pressure sources all around it.) In practice, a suitable grid is imposed, shot and receiver locations are assigned to gridpoints, and reconstruction is performed for each gridpoint in the grid. The reconstructed gathers are stored on disk.        2. For each trace, densely sampled gathers corresponding to the shot and receiver position of this trace are read from disk, their co-located traces are convolved with each other, and the convolution results are summed (stacked) to produce a single trace.        3. The convolution result (a single trace) contains the desired multiple prediction and is saved to disk.        
FIG. 1A shows traces from two gathers, centered at the source (source positions are denoted by a star in the drawings) and receiver (denoted by a diamond), which are convolved at co-located bounce points (denoted by circles) and the results summed. These SRME implementations have been used successfully in 2D multiple prediction. 2D prediction typically considers only the data from a single CMP (common mid-point) line at a time. In contrast, 3D SRME requires simultaneous use of traces from many nearby CMP lines.
The main drawback of the method described above in the 3D case is the need to repeatedly read large amounts of data (each reconstructed shot gather can be as large as 1-3 Gb) from disk and the associated need to store all, or at least large portions of, the reconstructed data on disk (the total size of reconstructed data can reach hundreds of Tb for a typical 3D survey).
The difficulty of reconstructing and handling vast amounts of data involved has been recognized, and various ways of circumventing it have been proposed. Early attempts to improve efficiency of multiple prediction focused on parallelizing computation over temporal frequencies (see “Parallelisation of surface-related multiple elimination,” van Waveren and Godfrey, Springer Berlin/Heidelberg, High-Performance Computing and Networking 919, 156-163 (1995)). Parallelizing is a necessary feature of most modern implementations of convolutional data-driven multiple prediction, but it does not resolve the problem of having to store and repeatedly access large volumes of data. Other authors circumvent the problem by not explicitly reconstructing the data. See for example van Dedem and Verschuur, “3D surface multiple prediction using sparse inversion,” 71st Internat. Mtg. Soc. Expl. Geophys., Expanded Abstracts, 1285-1288 (2001). Their method attempts to recreate the missing data on-the-fly using certain assumptions about the character of the data. It requires higher than usual acquisition effort and does not work well with conventional streamer data. Moore et al. (PCT Patent Application Publication WO 2005/103764A1) also use on-the-fly data reconstruction by selecting the nearest (in terms of spatial location, offset, and azimuth) available field trace to compensate for missing data. This approach may prove to be inadequate in the presence of rapidly changing complex structure, where the nearest existing trace may not be a good approximation to the missing one. Matson et al. in U.S. Pat. No. 7,197,399 (“Method of multiple prediction”) propose to reconstruct 3D multiple predictions from 2D predictions instead of reconstructing the data necessary to obtain the 3D predictions directly. This method requires the use of substantially different velocities for multiples of different orders and may be insufficiently accurate when multiples of different orders overlap. The latter is often the case, especially in interbed multiple prediction.
An improved method is needed to deal with the difficulty of reconstructing and handling vast amounts of data in 3D SRME or IME. The present invention satisfies this need.